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SPECTRAL TRANSFORMATIONS IN THE "SOFI" COMPLEX FOR 
PROCESSING PHOTOGRAPHIC IMAGES ON THE ES COMPUTER. PART I 

A. S. Debabov and D. A. Usikov 


INTRODUCTION /3* 

The study [1] determined the tasks and pirnciples for 
organizing a system of mathematical processing of photographic 
images (the "SOFI" complex). Spectral transformations play no 
small role in the programming support for the complex. Since 
the information is inputted to the computer in discrete form, 
we shall be speaking of discrete transformations (DT). 

The DT of ;he numerical file C(NNN) represents a series of 
linear transformations of the subfiles f^ N \ selected from the 
f ile according to a definite rule: 


u 



* m 0, /,... A/- /; Af##, 


( 1 ) 


where K is an integral function of k; N — DT base. Drawing an 
analogy between the subfile and the vector column in N-dimen- 
sional space, we may write: 
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to .I'/r'j** f wt 

•F £ a . t , 
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( 2 ) 


Mark ° over the sympol of the subfile means that its °-tranS" 
forms, i.e. linear combinations of the subfile elements, occupy 
the position of the subfile elements in the computer memory. 
f| and q — symbols of the DT matrix of the dimensionality N 

Numbers in margin indicate pagination of original foreign text. 
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x N and its elements; p — normed factor. In the general case 
to implement DT (2), we require ~ N** operations of addition- 

a 

multiplication (’ **) , and a memory for - tr numbers — elem- 

ents of the DT matrix. 

Since 1965 , methods have been vigorously developed for the /Jj_ 
rapid transformations (RT) [2, 3,^,5], based on the possibility 
of representing the DT matrices IT ^ in the form of Kronecker 
products of the slightly filled matrices ■ L,2,...,m). 

In this case, (2) is transformed into the DT series 


where 


h m r 


f 


«*» 


f*Lf» 

r * *M 



(3) 


if cf**- 4 - - are block-diagonal matrices, the number of 

operations £* 1 ) and the volume of additional memory are 
reduced by a factor of t v \ , where 



"Acceleration" of the DT on a computer Is achieved by: 

x - ("• f ) - (*. m J 

a) finding ‘■he factorisation method H ~Q ...0 

for which i (formula (4)) is maximum; 
rt 

b) constructing the optimum algorithm for the DT series 
(formula (3)) with allowance for the layout of the subfile in 
the computer memory (formula (1)); 

c) optimum translation of the algorithm into a system of 
actual computer commands; 
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d) technical improvements in the computer: special 

commands, special processors, etc. 

Algorithms and programs for one-dimensional Poui-ier RT 
(RTF) are the most widely used at the present time. The mod- 
ule of the multidimensional RTF described below is ”mo.*e rapid" 
than its predecessors In all points; with respect to (d) we 
have in mind the advantages of the EC computer as a third 
generation macnine; with respect to a), b) — the advantages 
of the RTF are a real-valued function [6], extraction of a 
normed multiplier, inversion, equal multipliers, exception of 
transposition (A. S. Debabov), the extraction of 0.1 multi- 
pliers, and recurrent computations of trigonometric functions 
(N. I. Kalinin); with respect to c) — programming on ASSEM- /5 
BLER of cyclical DT procedures, the selection of optimum 
command sequences based on results of tests for rapid operation 
(A. S. Debabov). The DTF of a file of ** 10^ numbers on one or 
several measurements on the EC-1040 computer is performed in a 
time of approximately 10 seconds. 

For several applications, the Fourier transforms of a real- 
valued file should be presented in the form of a complex file, 
because from the latter we must extract the purely real and the 
purely imaginary transforms, etc. For a two-dimensional file, 
a special subprogram performs similar transformations, which is 
described below. The third subprogram makes it possible to 
calculate the image correlation functions with the help of the 
first two. The example given here may be regarded as the text 
of the module described. 


I. STRUCTURE OF INFORMATION IN THE COMPUTER MEMORY 

In accordance with the principles for organizing the 
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"SOPI” complex [1], the exchange of information between the 
user programs and the DT modules is performed by using the 
unlabelled general region of the form 


COMMON nCT(J9*),NC,N,##f,NA^*, C(MNQ... , ^ 0 ) 

where FICT is the region of the control program [1]; «NC 
the control parameters of the DT modules described and the 
auxiliary routines; NNN — total dimensions of the file of 
numbers of the HEAL * 4 type. 

During the operation of the user program phase with the 
DT modules, requiring a memory of Z* t bytes, and Zp * bytes, 

HNN < 256 ( Zp - Z fl ) Thus, if Zp m &** , and * *** , 

the user may write up to 15 files of the form A2(64,6 i O, or 
up to 3 files of B2( 128,128), or a three-dimensional file C3 
(32,32,32), etc. In the operator (5) producing the program in the memory. /6 
Some files may be first written as complex lengths of 8, for example, CMP 
(128^128), which is equivalent to the real file description 03(2,128,128). 

In the DT subroutine, the file C(NNN) may be written with 
its first element C ( I ) and the rules (I) for selecting the 
subfiles, which are determined by the values of the parameters 
NC-NK. As an example, we give the rule for selecting the sub- 
files on three measurements of the file C3(l*2,3), which are 
in operation in the FPT modulus at INC I * I : the parameters 
NA and NK determine the boundaries of C3 in the file field C(NNN): 
C3(I , I , I , )=C (NA ) , C3 (NI ,N2,N3 )=C (NK) ; the parameter N determines 
the base of the first measurement (row of two-dimensional 
subfiles): NI - N, with the selection of the k-thelement of the 
i-th subfile of the first measurement (0 s < K <N), (0 N2 x N3) 

(I) has the form: K = WA+k + i*Af ; the parameter NN determines 

the base of the two-dimensional subfiles first and second 


k 


measurements: N2 - NN/N, (1) for the k-th element of the j-th 

subfile of the second measurement iOsk^/Zc), (0 *j<N ) , selected 

from [1] of the a— t h two-dimensional subfile such that: 

K*^A+ , where os( < yVJ ; the base of the third 

measurement (numbers of the two-dimensional subfiles N3 ■ 

NXtl-NA )/NN » anci the ru le (1): , where 

Qfm <NN 

The transforms of the subfiles of the corresponding measure- 
ments are selected according to the same rules as the subfiles 
until DT is performed. 

In the INRECO module, the elements of the real two-dimen- 
sional subfile C2 are sele ted according to the rules; 

C 2 U,;)*Cl/VAti , and the corresponding elements of the com- 
plex subfile t), f where 

u$l<N% oSjk corresponds to the real part; e ■ 1 — 

the imaginary part. 


2. RTF ALGORITHM OF A REAL-VALUED FUNCTION /7 


The canonical form of the DTF record of the complex 
function r° of the discrete whole-number argument k on the 
base (0 #KAM) [2] is: 

<«• (6 

and thus the inverse transformation: 


r-c<;r| 


M "l". 


t 


&cp<tr 


( 7 ) 
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where ) • 00$ t $n )• ^ N m < c ti • >• 


the value of the upper Index for W N 1 b the exponent, and for 
and — the multiplicity factor of the angle to the base 
angle ^ , oj • 00 » ( ^ >. 


The multidimensional DTF may be presented in the form of a 
sequence of one-dimensional DTF; for the function of the two 
arguments on the base a/A/«A/1+^2i 



(80 


since 

Substituting Im. r-o into (7) and setting Re{£Ur' 
for odd N, we obtain (and we set): 


Re fT'-Refe-# , . 

( ( 9 ) 


In this notation, the conversion formula takes the form: 



(10) 


i.e., the cosine and sine coefficients of the Fourier series /8 
may be written toward to each other. 

We may obtain the algorithm of the rapid reverse DTF for 
M 

N = 2 by comparing the expansion of (10) for odd and even 
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elements of the subfile with the expansion (10) of 

auxiliary subfiles and , selected according to 


the rule: 


For example, for even elements for all h -It 


Considering that 
we obtain 


< •< . S2 . SI . e|*'.ejt 


I. 

t. 

3 . 

4 . 


* 


ir-o 

S 0 ♦! 


<*» £<»►) 

i L*?i . 8 ri * 


« 


b *-j oi,-j “ * 


Performing similar, although more cumber- 
some caleuations for the odd elements, we obtain: 


where 


tr-e*-«r, 

* *• > 

tf-ktizwi), 

jw-ir-fc; . 


s. 

4 . 

7 . 

6. 


(12) 


K -It 


Formulas (12) hold for any h.*i , and at h * 1: 


/9 


I?. ft. . 
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te* fe i 

Assuming *1 ••» for n ■ N/2 and applying (12) in 

the case to sequentially ordered auxiliary 

subfiles, whose number doubles with each step, at the locations 
of the elements jf** we obtain the shuttled elements of the 

subfile d$> , i.e., the DT result (10). The number of the 

element is connected with the number of the element k « l 
by the law of binary inversion: 

if 


m * 0 * 

♦ ij2” ♦ • •• ♦ * 11 - 2 * ♦ ( 13 ) 

then 

where i-I) assumes values of either 0 or 1. For 

brevity, we shall call the procedure and the result of trans- 
posing the subfile elements according to the rule (13) as 
inversion . 


Direct RTF of tht real-valued function: 



(14) 


produces the DT sequence described above (including inversion), 
accomplished in reverse order. The algorithm for changing 
from the subfiles >> t* to the subfile is pro- 

duced by solving a linear system of equations (12) relative 

to If . 

The number of arithmetic operations for RTF(10) at M>3* 
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and f and with the Indirect use 

of the algorithm (10) — [*(*'«)]•*’ and */*♦«)] ft* . 

During transformation of the m-dlmenslonal files, having the 
Identical base N In all the measurements , the number of opera- 
tions Increases by a factor of ( m ■ 3 » 3 . . . ) • Direct 

RTF (Hi) requires as many as f*« and [N-2] fewer «+9 than 

the reverse RTF (mulitpllcation by two, we refer to operations 
of the type of ' ♦ *), If we let the user combine multiplica- 
tion of each transform element by I/N (In the multidimensional 
case — by i/( 4 / 1 *A/2+... */im )) with the usually necessary 

multiplication oy values calculated previously, which we do. 

If the user is Interested In complex representations of 
the transform (6), In the one-dlmenslonal case, we may change 
to It by using the formulas (9). For man> measurements, the 
change algorithm Is complicated. Thus, the two-dimensional DTF 
of the real-valued function OC* 

1 e r, r(«) 

. (15) 

Comparing (15) and (8), at Cl -0 ,VV2 and at ll • 0 ,to /2 we 

A 

obtain the formulas for connecting ^ and n of the type 
(9), and at I < Cl « A/I/2-1 , U t2< Hl/l-U 


A 


(Ml, Ml) 

fi.Ci 



i^Ml.Ml) 



-f 


A»i-h,U * 



+ 






; - ; i, 

i }. 

1 1 

i ♦ ♦ ]. 


( 16 ) 


/10 
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The formulas (9)» ( 16 ), their inversions with respect to f” 

and analogous formulas for the case Rt P m 4), A ■ 1m F represent 
the algorithm for the subprogram INRECO described below. 

3 . PPTD SUuPROGRAM 

The PPTD module performs direct (14, 15) and reciprocal 
(10) DTF for the real-valued subfiles (REAL # 4 ) of one ir 

M 

several measurements with the bases m * 2 (M is a whole number) 
according to each measurement or accoruing to several measure- /II 
ments. It is assumed that the processed files are fully loca- 
ted In the computer operational memory . Their Fourier trans- 
forms (increased by a factor of N) occupy the locations of the 
subfiles in the case of direct transformation; reciprocal trans- 
formation restores the subfiles (within an accuracy of the 
factor N) . 

The DTP is performed by using the RTF algorithm, described 
in the previous section (12, 13). Inversion is performed by 
separate transformation to PPTD at NC ■ 0. The rules for 
selecting the subfiles were described in Section I. The sub- 
files selected according to the general law were developed in 
parallel. The cases J ati/to ; flj* ■ were extracted 

in the algorithm (12), which leads to a reduction in the 
number of operations. 

The text of the PFTD subprogram in the FORTRAN-4 language 
is shown in Appendix 1. The PFTD module, cataloged in the 
library of target modules of the "SOFI" complex [1] in the 
Institute of Space Research, Academy of Science USSR retains 
the logic of this text, but the binary cycles 11, 22, 29 are 
replaced by invocations of the PPM module, which is written 
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In ASSEMBLER with allowance for the results of the tests on 
high-speeds groups of machl.ie commands performed on the ES- 
1040* One result of this study was an Increase in the actual 
operation of the ?-TD by a factor of 3.5» as compared with 
the FORTRAN option. As an example , the direct and reciprocal 
two-dimensional RTF of the file B2(128,128) require the 
together of 5-6 seconds ES-1040 processor time whereas the 
translation of the text (Appendix 1) on the dESM-6 computer 
requires 20-30 seconds to perform these procedures. 

At M,s 10, the normed result of the direct reciprocal DTF 
coincides with the initial (in terms of all the elements) 
subfile (which differ by 5 orders of magnitude in terms of the 
module), with a relative error of < 10“^. When these pro- 
cedures are repeated many times, the error increases linearly. 

At M > 10 the error increases due to the Inexact values of 
the factors °L > sL , which are not stored in the memory 
in this variation of the RTF, but are calculated recurrently 
from 03 * t in the DT process. /12 

The module call from the user program, written in FORTRAN, 
is performed as follows: CALL FFTD. Before the call, it is 

necessary to write the total region by the operator of the 
type (5), assign the necessary values of the parameters NC, 

N, NN, NA, NK and fill the processing section of the file 
C(NNN) with the necessary information. 

FFTD takes into account the values of the five parameters 
INTEGER * 4: NA and NK -- respectively the number of the 

first and last elements of the processed group of subfiles in 
the file field C(NNN) : U A/A<WK*A/AW i 


N — base of the subfiles of the first measurement: N1 - N : 
NN — base of two-dimensional subfiles ol the first and second 
measurements: N2 - NN/N; 
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NC — splitting parameter: at NC ■ 0 only the Inversion of 

subfiles 0 . first measurement is performed, if NN/N < 16 ; in 
the oppoelte case — inversion of subfiles of the first and 
second measurements and if (/* *V * 16 , the section is 

subjected to inversion for three measurements: NC > 0 
direct RTF inversion; NC < 0 — reciprocal RTF: 
at 

|^CJ ■ 1 - using three measurements 

|A/C| - 2 - using the first two measurements, 

Ja/CI ■ 3 - using the first measurement (rows) 

J fJC I ■ - using the second measurement (columns), 

|A/C| ■ 5 - using the third measurement (number of pages). 

The memory used by the module is 2.8 K. 

4 . INRECO SUBPROGRAM 

The INRECO module accomplishes the exchange of data between 
real - A2 and complex - CMP two-dimensional files with identi- 
cal dimensions, determined by the user on ove:'apping segments 
of the file C(NNN). In one of the files, eit! *• the function 
of two discrete arguments or its transform (8) or (15) is /13 

given . 

Access to the FORTRAN program: CALL INRECO. The para- 
meters: NC - NK. The meaning of the parameters N, NN, NC is 

the same as for FFTD (see Sec. 3); NA -- number of the first 
element of the r°al file: NK — number of the real part of 
the first element of the complex file in the C(NNN) field; 
the splitting parameter NC assumes values from 1 to 8, for 0 * 
dd values of NC XX*FE — procedures with the real part of the 
complex function; for even NC XX = IM — procedures with the 
purely imaginary part: 

at NC * 1 the XX transform is carried from A2 into CMP; 
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at NC ■ 2 the XX transform from A2 is added to the CMP: 

NC ■ 3,4 the XX transforms are extracted from CMP and 

moved to A2; 

NC ■ 5*6 the A2 1 b carried to the XX part of CMP; 

NC ■ 7*8 the XX part of CMP is carried to A2. 

Thus, operations at NC < 4 pertain to the Fourier image spectra; 

at NC 5 — to the images thei jives. 


The occupable memory is 3.0 K. 


5 . FFT2 SUBPROGRAM 

The FFT2 module combines three frequently encountered DT 
procedures of two- dimensional files, determined in the user pro- 
gram as complex: 

1) direct DTF of the two-dimensional file (image) with a 
zero imaginary part; 

2) inverse DTF of the complex two-dimensional Fourier 
transform (it was not apparent previously that the 
imaginary portion of the result is null); 

3) multiplication of the complex image spectrum by the 

complex conjugate of the "reference" spectrum. /l4 


Each procedure is carried out for individual access to 

FFT2 : 

1 - at NC > 0, 

2 - at NC < 0, 

3 - at NC - 0. 

For procedures 1 and 2, the user must determine on the file 
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field C(NNN) one or two real filer of the same dimensions ar the 
complex ones: NI ■ N, N2 ■ NN/N; NK designates the origin of 

the complex file; and NA designates the position of the first 
element of the first of the working real files. For procedure 
3, the user determines two complex files of the dimensions given 
above, whose initial address corresponds to the C(NK) and C(NA) 
addresses; the conjugate is given according to elements of the 
second, and the result is written in the first file. The second 
complex file and the working real files may be determined in the 
same section of the field (CNNN) . 

FFT2 calls the FFTD and INRECO modules. It may 3erve as an 
example of access to this sub. The FFT2 text is presented in 
Appendix II. The occupa'ole memory is 0.9 K. 


6. COMPUTATION OF CORRELATION FUNCTIONS. EXAMPLE 

An example of using the moduli described above, and also 
their text, is the program whose text is presented in Appendix 
III. The bases of e files (16, 16 ) were selected for ease 
of visualization of the image function and the correlation 
function, performed by the PR2 module. The KREST module per- 
forms generation of the reference and the images. Print-outs 
from the PR2 are given in Appendix IV: 

a) The reference-cross with the center coordinates (6,8). /15 

The coordinates are selected to the right and below che ele- 
ment (1,1). In view c-f the periodic nature of the images, 
for negative values of a certain coordinate the origin of the 
coordinate will be one of the positions designated by a cross; 


b and c) Images No. 1 and NO. 2 represent two crosses with 


the parameters of the reference, shifted with respect to the 
last and with respect to each other. The dotted line combines 
the cross contours. The arrows are vectors of the reference 
center to the center of the image element; 

d and e) The correlation functions of the images No. 1 
and No. 2 with the reference a), obtained as the result of 
calculations using the program (Appendix IV). The arrows are 
the vectors from the "corresponding" origin of the coordinates 
to the position of the correlation function maximum; the dashed 
line is drawn on the level of the values of the correlation 
function module, which equals half of the maximum value. 

In conclusion, the authors would like t <-> thank N. I. 
Kalinin, T.E. Krenkelya and B. I. Kolosov for valuable dis- 
cussions . 
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